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ABSTRACT 

We investigate the evolution of substructure in cold dark matter halos using N-body simulations 
of tidal stripping of substructure halos within a static host potential. We find that halos modeled 
following the Navarro, Frenk & White (NFW) mass profile lose mass continuously due to tides from the 
massive host, leading to the total disruption of satellite halos with small tidal radii. Although mass is 
predominantly stripped from the outer regions, tidal heating also causes the halo to expand and the 
central density to decrease after each pericentric passage, when mass loss preferentially occurs. As a 
result, simple models based on the tidal- limit approximation underestimate significantly the tidal mass 
loss over several orbits. The equilibrium structure of stripped NFW halos depends mainly on the fraction 
of mass lost, and can be expressed in terms of a simple correction to the original NFW profile. We apply 
these results to substructure in the Milky Way, and conclude that the dark matter halos surrounding 
its dwarf spheroidal (dSph) satellites have circular velocity curves that peak well beyond the luminous 
radius at velocities significantly higher than expected from the stellar line-of-sight velocity dispersion. 
Our modeling suggests that the true tidal radii of dSphs lie well beyond the putative tidal cutoff observed 
in the surface brightness profile, suggesting that the latter are not really tidal in origin but rather features 
in the light profile of limited dynamical relevance. For Draco, in particular, our modeling implies that 
its tidal radius is much larger than derived by Irwin & Hatzidimitriou (1995), lending support to the 
interpretation of recent Sloan survey data by Odenkirchen et al. (2001). Similarly, our model suggests 
that Carina's halo has a peak circular velocity of order ~ 55 kms -1 , which may help explain how this 
small galaxy has managed to retain enough gas to undergo several bursts of star formation. Our results 
imply a close correspondence between the most massive substructure halos expected in a CDM universe 
and the known satellites of the Milky Way, and suggest that only substructure halos with peak circular 
velocities below 35 km s" 1 lack readily detectable luminous counterparts. 

Subject headings: N-body simulations; galaxies: evolution, structure; cosmology: dark matter 
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1. INTRODUCTION 

Over the past two decades, cosmological models based 
on the assumption that the Universe is gravitationally 
dominated by Cold Dark Matter (CDM) have evolved 
into the prevailing paradigm for interpreting the forma- 
tion and evolution of structure in the Universe (Blumen- 
thal et al. 1984, Davis et al. 1985, Bahcall et al. 1999). In 
this model, structure develops as primordial density fluc- 
tuations with scale-invariant power spectra grow to form 
virialized structures. Due to the negligible thermal ve- 
locity of CDM, fluctuations survive the early universe on 
all scales, and structure develops from the bottom up as 
small, dense CDM clumps collapse at very early times and 
subsequently undergo a series of mergers that result in the 
hierarchical formation of large, massive dark matter halos. 

These halos constitute the hosts of all galaxy systems, 
from individual galaxies to galaxy groups and clusters. In 
a pioneering paper, White & Rees (1978) noted that merg- 
ers of dark matter clumps should very efficiently disrupt 
the progenitors and lead to smooth remnants where most 
of the mass is in a single, monolithic structure. Early cos- 
mological N-body simulations showed that the remnants of 
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the merger hierarchy were featureless systems with little 
substructure, and appeared to confirm this expectation. 
The survival of individual galaxies in groups and clusters 
was thus linked to the ability of baryons to cool and con- 
dense into much more tightly bound structures than their 
surrounding dark halos, allowing them to survive as self- 
bound entities the assembly of galaxy clusters. 

Recently, the development of efficient algorithms and 
the advent of massively-parallel computers have led to sim- 
ulations with dramatically improved force and mass res- 
olution. These have revealed a wealth of substructure in 
virialized dark matter halos, where none was apparent at 
lower resolution (Ghigna et al. 1998, Klypin et al. 1999b, 
hereafter K99b, Springcl et al. 2000). Substructure ap- 
pears to extend down to the smallest mass scales reliably 
resolved, and several hundred dark matter satellites are ex- 
pected within a few hundred kpc of the Milky Way. Most 
of these substructure halos, hereafter 'subhalos', are of low 
mass. Indeed, as a whole, subhalos contribute only about 
10% of the mass of a virialized halo, and most of the mass is 
in a smooth monolithic structure, as envisioned by White 
& Rees (1978). 
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These results imply that CDM halos are quite resilient 
to tidal disruption, since they may lose a large fraction of 
their mass without being completely disrupted. This re- 
silience stems from the large concentration of their mass 
distribution, which is in turn intimately linked to the neg- 
ligible thermal speed of CDM particles. Indeed, models 
where dark matter is 'hot' or 'warm' rather than cold lead 
to the formation of halos with much reduced substructure 
(Colin et al. 2000, Bode, Ostriker & Turok 2001, Avila- 
Reese et al. 2001, Dalcanton & Hogan 2001, Knebe et 
al. 2002). 

The dependence of substructure on the nature of dark 
matter makes it a discriminating tool between competing 
cosmogonical models. For example, Klypin et al. (1999a, 
hereafter K99a) and Moore et al. (1999, hereafter M99) 
have pointed out that the number of low mass subhalos 
identified in CDM simulations of a Milky Way-sized galaxy 
halo exceeds the number of observed Galactic satellites by 
a factor of 10 to 100. This observation has fueled spec- 
ulation that it may be necessary to modify radically the 
CDM paradigm on small scales in order to account for the 
scarcity of luminous satellites in the vicinity of the Milky 
Way. 

Possible modifications include a finite interaction cross- 
section for dark matter particles (Spergel & Stcinhardt 
2000), or a finite 'temperature' for dark matter, as in the 
warm dark matter scenario (Dalcanton & Hogan 2001). 
Whether or not such radical modifications are actually 
necessary is still being debated, although it remains pos- 
sible that the apparent 'decoupling' between dark and lu- 
minous matter in low-mass halos may be a natural con- 
sequence of the effect of stellar feedback on these scales 
(White & Rees 1978, Kauffmann, White & Guiderdoni 
1993, Cole et al 1994), or perhaps of the heating effects 
of an ionizing UV background (Efstathiou 1992, Bullock, 
Kravtsov & Weinberg 2000, Somerville 2001, Benson et 
al. 2001). 

So far, most of the effort devoted to studying the prop- 
erties and evolution of substructure have focused on a few 
cosmological simulations with ultra-high numerical reso- 
lution, such as those performed by Moore et al. (1998), 
Ghigna et al. (2000) and Klypin et al. (2001). Such stud- 
ies are hindered by the small fraction (<~ 10%) of mass 
that remains attached to self-bound subhalos in a typical 
virialized system. As a result, even in the best simulations 
a typical subhalo is made up of little more than ~ 10 3 par- 
ticles and is particularly vulnerable to numerical artifacts 
affecting small-N systems near the resolution limit. In ad- 
dition, the complexity of the cosmological context, with its 
continuous mergers and rapidly changing potential wells, 
makes it difficult to isolate the mechanisms responsible for 
driving the evolution of substructure. Many important 
questions thus remain unanswered: for example, (i) are 
substructure halos fully disrupted by the tidal field of the 
host system?; (ii) how does the internal structure of sub- 
structure halos evolve as they are tidally stripped?; (iii) 
how do numerical limitations affect the mass function of 
subhalos? 

Addressing these questions requires much better res- 
olution than is achievable in current cosmological sim- 
ulations. We have therefore decided to investigate the 
evolution of subhalos by performing high-resolution sim- 
ulations of the evolution of realistic dark halo models 



within the gravitational potential of a much larger system. 
The model we adopt assumes a spherically-symmetric and 
time-independent host potential, but accounts fully for 
the self-interaction between particles in the subhalo. Our 
models are thus useful to study in detail the process of 
tidal heating and stripping by the main potential but ne- 
glect the complicating effects of dynamical friction, of ha- 
rassment between subhalos, and of severe departures from 
spherical symmetry in the host. Despite these caveats, our 
simulations can be usefully applied to the evolution of sys- 
tems much less massive than the host, as is the case for 
individual galaxies in a rich cluster or for satellites orbiting 
a bright galaxy. 

The outline of this paper is as follows. In Section 2 
we present our halo model and discuss the properties of 
equilibrium N-body models evolved in isolation. Section 3 
describes our simulations of satellite halos on bound orbits 
in the gravitational potential of a static host. Section 4 
compares the results of these simulations with the predic- 
tions of simple models based on the tidal limit and on the 
impulse approximation. Section 5 examines the changes in 
the internal structure of tidally stripped halos and Section 
6 explores its implications for the properties of dark halos 
surrounding the dwarf spheroidal companions of the Milky 
Way. We summarize our main conclusions in Section 7. 

2. THE HALO MODEL 

The survival of a dark matter halo in the presence of 
strong tides depends sensitively on its internal structure. 
It is therefore important to choose a halo model that repre- 
sents as faithfully as possible the structure of dark matter 
clumps expected to form in the cosmological model to be 
investigated. The density profile of cold dark matter halos 
has been the subject of extensive numerical studies by nu- 
merous authors (Frenk et al. 1988, Dubinski & Carlberg 
1991, Warren et al. 1992, Navarro, Frenk & White 1996, 
1997, Fukushige & Makino 1997, Moore et al. 1998), all 
of which agree on two basic results: (i) the density profile 
of CDM halos is shallower than isothermal near the center 
and steepens gradually outwards, becoming steeper than 
isothermal in the outer regions; and (ii) there is no well 
defined value for the central density of the dark matter, 
which can in principle climb to arbitrarily large values near 
the center. 

Navarro, Frenk & White (1996, 1997, hereafter NFW) 
propose a simple formula that synthesizes this result 
and that provides a good description of the spherically- 
averaged density profile of virialized CDM halos, 

PNFwM = —, , PS . rj, (1) 

r/r s (l + r/r s y 

where r s is a characteristic scale radius. Although there 
is some disagreement regarding how accurately this sim- 
ple formula describes the latest high-resolution simulations 
(Fukushige & Makino 1997, Moore et al. 1998, Ghigna et 
al. 2000, Klypin et al. 2001), the enclosed mass profile 
it predicts seems to agree well with high-resolution sim- 
ulations of galaxy-sized systems (Navarro 2001, Power et 
al. 2002). In view of this, and given its simplicity, we adopt 
this model for our study. 

One important feature of models such as NFW is that, 
when truncated, the total energy of the remaining sys- 
tem is not necessarily negative. This is shown in Fig- 
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ure 1, where we plot the kinetic and potential energy of 
an isotropic NFW halo truncated 'instantaneously' at ra- 
dius r (solid lines). The total energy of the truncated sys- 
tem is negative only when the truncation radius is larger 
than a 'minimum binding radius', r"bmd — 0.77 r s . This 
suggests that NFW models tidally limited to radii compa- 
rable to rbind (or, equivalently, that are stripped of more 
than ~ 91% of their mass) may be totally dispersed by 
tides. 

The dashed lines in Figure 1 show that this is not the 
case for a truncated singular isothermal sphere, which re- 
mains self-bound regardless of the truncation radius. As 
discussed by Moore, Katz, & Lake (1996), halos with sin- 
gular isothermal density profiles "should always survive at 
some level," but halos with shallower inner density profiles 
like NFW may disrupt unless a dissipational component is 
included. Our simulations are designed to study the mass 
loss timescale for halos orbiting within more massive sys- 
tems, as well as the conditions under which disruption by 
tides may occur. 

2.1. N-body Realizations of NFW Halos 

Our procedure for initializing isotropic N-body realiza- 
tions of NFW models follows closely the prescriptions of 
Hcrnquist (1993). Since NFW models have infinite mass, 
it is necessary to truncate them at some fiducial radius, 
which we take to be r cut = 10 r s . With this choice, 
the half-mass radius of the model is th ~ 3.6 r s . We 
choose the spline softening length of our runs, e, according 
to the scaling suggested by van Kampen (2000), so that 
e = 0.77 rh N~ x ^, where N is the total number of par- 
ticles. With this scaling, e = 0.06 r s and e = 0.19 r s for 
realizations with N = 10 5 and N = 3000, respectively. 3 

Unless otherwise specified, we shall adopt units where 
the gravitational constant, G, the scale radius, r s , and the 
total mass of the satellite halo, m sat , are all unity. With 
this choice, the crossing time of the system is t croS s = 
V r cut/^ m sat = 31.6 and the circular orbit timescale at 
the half-mass radius is t c irc(f/i) = 27rr^/V' c (r^) = 60.7. All 
simulations have been carried out using Stadel & Quinn's 
multi-stepping, parallel treecode PKDGRAV (Stadel 2001). 
We use a kick-drift-kick leapfrog integrator with adap- 
tive, individual timesteps chosen according to the gravita- 
tional potential and the acceleration, At, = 0.03 v^i/ffli- 
Timesteps in isolated runs typically range between 2.5 x 



10 4 and 4x10 3 i cross , with a median value of 10 
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mum roughly at r s , and decreases further out. This be- 
haviour is mimicked by the circular velocity, which peaks 
at r max = 2.16 r s . This radial dependence of the 'temper- 
ature' of the system determines to a large extent the main 
qualitative features of the relaxation-driven evolution of 
an N-body NFW model. 

The solid (dotted) curves in Figure 2 show the mass 
within various radii for an N — 3000 (N = 10 5 ) model 
evolved for up to 400 (80) crossing times. Two differ- 
ent evolutionary stages can be seen over this very long 
timescale. At first (t <J 100t croS s), the region within one 
scale radius (which is originally cooler than its surround- 
ings) is heated through collisions with faster-moving par- 
ticles and gradually expands to accommodate the energy 
gained. This is shown in Figure 2, where the mass enclosed 
within radii smaller than r s is seen to decrease systemati- 
cally for t < 100 icross- A similar process appears to oper- 
ate well beyond r s , which is seen to gain energy (i.e., to be- 
come less dense) at the expense of the hotter, inner regions, 
albeit on a much longer timescale (see, e.g., the curve la- 
belled 8.0 r s in Figure 2). To conserve energy, intermedi- 
ate regions of the system become increasingly bound and 
dense: the mass within 2 r s is seen to increase steadily 
throughout the evolution. We have explicitly checked that 
the evolution depicted in Figure 2 is not affected by inte- 
gration errors; our conservative timestepping choice leads 
to energy conservation better than 0.32% even after 400 
crossing times. We have also ruled out non-equilibrium 
initial conditions as the cause of the observed evolution. 
After some initial fluctuations within the first few crossing 
times, the ratio of kinetic to potential energy of the halo 
remains constant to within less than 2% for the duration 
of the simulation. 

For N = 3000, energy transfer into the inner core shuts 
off at t w 100 t cr oss, when the velocity dispersion in the 
inner regions ('core') has been raised to roughly match 
that at r s . The second stage of the process begins then, 
in which the core contracts as more and more particles 
are kicked into highly energetic orbits. This 'core collapse' 
stage, familiar from globular cluster evolution calculations, 
leads gradually to highly concentrated systems: after 400 
crossing times the density within 0.25 r s has gone up by 
almost a factor of 3 and the mass within r s has nearly 
doubled. 

On what timescale does the inner mass profile evolve? 
The bottom panel in Figure 3 shows the mass within 0.025, 
0.05, 0.1, and 0.2 r s , respectively, for a run with N = 10 6 
particles. The mass within these radii (shown in units of 
the theoretical mass within each radius No(r)) decreases 
continuously and at rates that increase gradually nearer 
the center. We measure the timescale of this evolution 
from fits of the form, N(t) — N a e~*/*°, shown as straight 
solid lines in Figure 3. These timescales are listed in Ta- 
ble 1, and exceed by about two orders of magnitude the 
local collisional relaxation timescale, 



2.2. Secular Evolution of N-body Halo Models 

Isolated N-body halo models evolve away from their 
equilibrium configuration due to discreteness effects as- 
sociated with the finite number of particles. The main 
mechanism driving this secular evolution is encounters be- 
tween particles, which lead to energy exchanges between 
different parts of the system and, consequently, to gradual 
modifications in the structure of the halo. 

It is instructive to consider the specific kinetic energy 
of particles at various radii throughout the system, since 
collisions will naturally tend to drive the system towards 
cquipartition. The velocity dispersion in isotropic NFW 

halos rises from the center outwards, reaches a maxi- Indeed, Table 1 shows that mass is evacuated from the 

3 Dehnen (2001) suggests an alternative scaling for the optimal softening length, e oc N 1 / 5 for compact softening kernels. With this scaling, 
e = 0.12 r s and e = 0.21 r s for halo realizations with N = 10 s and N = 3000, respectively. 
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center at rates consistent with the evaporation timescale, 

icvapO) ~ 136i ro la X (r) (3) 

(Binney & Tremaine 1987). This is not entirely unex- 
pected; individual particle energies do change on the relax- 
ation timescale, but it is the slower process of evaporation 
that leads to appreciable systematic changes in the mass 
profile of halo models. The straight lines in the top panel 
of Figure 3 confirm that the same considerations also de- 
scribe quite well the evolution of the inner mass profile in 
the 3000-particle runs: these are not fits to the data, but 
rather show the decline in mass expected from the local 
evaporation timescale computed from eq. 3. 

We note that the secular evolution of spherically- 
symmetric N-body systems has been investigated previ- 
ously by Quinlan (1996) using Fokkcr-Planck methods. 
Quinlan's study indicates that collisional effects start mod- 
ifying the mass profile on the local relaxation timescale, 
whereas ours suggest that, once started, the mass profile 
gets modified at a rate controlled by the local evaporation 
timescale. Quantitatively, our results are in reasonably 
good agreement with Quinlan's work. For example, scal- 
ing his results to our NFW models with 10 5 particles, we 
would expect the mean density within a radius containing 
roughly 1% of the total mass to decrease by about 50% 
over a timescale equivalent to 100i cross . This is in rea- 
sonable agreement with the results presented in Figure 2, 
where the mean density within a fixed radius r = 0.25 r s 
(initially containing 1.6% of the total mass), is seen to 
decrease by about 30% over a similar timescale. 

In summary, we expect N-body models of NFW halos 
to deviate significantly from their original structure on 
the local evaporation timescale. Halos tend to develop an 
isothermal "core" of size r ~ r s on a timescale comparable 
to t eva , p (r s ). Over longer periods, the central regions are 
expected to become significantly denser, and thus more 
resilient to tides. This affects primarily small-iV systems, 
where relaxation timescales are short. 

We can estimate the importance of this process in cos- 
mological N-body simulations by identifying the outer ra- 
dius, r cut , with the 'virial radius', T2oo, where the mean 
density contrast 4 is 200. With this choice, t c i IC (r cut ) is 
comparable to the age of the universe at the time a halo is 
identified. Thus, halos formed at z = 10, 5, and 1, would 
have evolved (if kept isolated) for 570, 218, and 26 crossing 
times by z — 0, respectively. (We have assumed a ACDM 
cosmogony with £1 = 0.3, Oa = 0.7, h = 0.65). Accord- 
ing to eq. 3, the age of such systems at z — will match 
their half-mass evaporation timescales for N = 200,55, 
and N < 1, respectively. We conclude that it is unlikely 
that relaxation effects will introduce a significant bias in 
the evolution of substructure halos with N > 200. The 
survival of halos with A" <J 200 collapsing at z < 10, on 
the other hand, should be carefully monitored to rule out 
the possibility that they may be affected by relaxation ef- 
fects. 

3. TIDAL STRIPPING OF SUBHALOS 
3.1. The Numerical Setup 

4 We use the term 'density contrast' to denote densities expressed in 
the present value of Hubble's constant as H(z = 0) = Hq = 100 h km s 



We use the NFW halo models described in the previous 
section to investigate quantitatively the tidal stripping of 
'satellite' halos orbiting within a much more massive host 
system. To this effect, after allowing an NFW model to 
relax for 10 crossing times, we place it in orbit within a 
static, spherically-symmetric NFW potential with scale ra- 
dius R s = 10 r s , and total mass, Mhost, within 10i? s . Our 
simulations thus neglect the effects of dynamical friction 
and the reaction of the host to the satellite, but given the 
vast difference in mass and size of the two systems these 
are expected to contribute only minor corrections to our 
results. For example, applying the approximate formula 
of Colpi, Mayer & Governato (1999) we find that the dy- 
namical friction timescale for an orbit with r ap = 10i? s 
and r pcr = 1.5 i? s exceeds the orbital timescale by up to a 
factor of <~ 20. 

Note that due to the scale-free nature of gravity, the 
physical scaling of our systems is arbitrary. For example, 
if we take r s w 10 kpc and m sa t ~ 10 12 M Q the satellite 
halo would represent a Milky Way-sized galaxy, and the 
host potential would corresponds to that of a rich galaxy 
cluster. Choosing r s w 1 kpc and m sat w 1O 9 M , on the 
other hand, would correspond to the case of a dwarf galaxy 
orbiting within the potential of a Milky Way-sized halo. 

We have performed a series of simulations varying the 
orbit, the mass resolution of the satellite, and the mass of 
the host relative to that of the satellite. In particular, we 
have followed the evolution of A^ = 3000 and A" = 10 5 ha- 
los on three different eccentric orbits, all of which have the 
same apocentric radius, r ap = 10i? s , but different percen- 
ters, r pcr = 1.5, 3.0, and 6.0 R s , respectively. Note that the 
median apocenter to pericenter ratio of subhalos in cosmo- 
logical N-body simulations is between 5:1 and 6:1 (Ghigna 
et al. 1998, Tormen, Diaferio, & Syer 1998), comparable 
to some of our eccentric orbits (see Table 2). Several 3000- 
particle satellites were also placed on circular orbits in this 
potential, with orbital radii ranging from r c ; rc = 1.5 to 
30.0 R s . The mass of the host system is set to 300 m sat for 
these runs. 

Finally, we have investigated the conditions necessary 
for total disruption of an NFW halo by performing simu- 
lations oi N — 10 5 satellites on circular orbits with orbital 
radius r c i rc = 1.5 R s and varying the mass of the host sys- 
tem from Mhost = 300 m sat to M hos t = 3000 m sat . The 
numerical and orbital parameters of these simulations are 
listed in Table 2. 

3.2. Mass Loss 

Figure 4 shows a sequence of snapshots of an iV = 10 5 
satellite halo on an eccentric orbit with apocenter r ap = 
10 R s , and pericenter r pcr = 1.5 R s . The effects of tidal 
stripping and mass loss are clearly visible in the form of 
elongated tidal 'tails', which extend for many scale radii 
beyond the bound core of the satellite. The satellite ap- 
pears increasingly elongated along the path of its orbit 
as it approaches pericenter, and during each pericentric 
passage a substantial fraction of particles are physically 
separated from the satellite by the tidal field of the host. 

We quantify this mass loss by identifying the maxi- 
mal subset of particles which remain gravitationally self- 
bound. Our procedure follows that of Tormen et al. (1998), 

units of the critical density for closure, p cr it = 3H 2 /8nG. We express 
- 1 Mpc- 1 
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and identifies iteratively all particles that have negative 
binding energy relative to all others. At each iteration 
all escapers (i.e., particles with positive binding energy) 
are removed from the system and the potential and ki- 
netic energies of all remaining particles (in the rest frame 
of the bound subsystem) are recomputed. The iteration 
ends when no further escapers arc found. Note that this 
process is also used by the group finder SKID to iden- 
tify gravitationally bound groups in N-body simulations 
(Stadel 2001). 

The evolution of the self-bound mass, mb n d, is shown 
for three different eccentric orbits in Figure 5. The or- 
bital parameters are shown in each panel. The figure 
shows results for models with N = 3000 (dotted lines) 
and N = 10 5 (solid lines), and confirms our conclusion 
in § 2.2 that the mass loss process is rather insensitive to 
the mass resolution when the number of remaining bound 
particles exceeds a few hundred. Mass loss takes place 
in easily recognizable steps, where each event is triggered 
by tides operating at pericenter. As expected, the total 
amount of mass lost correlates strongly with pericentric 
distance; the satellite on the r pcr = 1.5 R s orbit is reduced 
to less than 5% of its original mass after 10 orbits, while 
that on the r per = 6.0 R s orbit retains 35% of its mass af- 
ter 10 orbits. Satellites appear to lose mass continuously, 
especially for orbits with small pericentric radii, suggest- 
ing that satellites are not in general stripped down to a 
finite tidal radius in a single orbit. Halos with the same 
pericentric distance may therefore be stripped to very dif- 
ferent degrees depending on the total number of orbits 
completed. 

How many orbits has a typical substructure halo com- 
pleted by z — 0? The orbital time depends mainly on 
the apocentric distance (see, e.g., Table 2), which may 
be assumed to scale with the time of accretion into the 
host system in roughly the same way that the turnaround 
radius scales with time in simple spherical infall models, 
r ta oc t 8 / 9 (Bcrtschingcr 1985). For example, assuming 
that, for systems accreted at z = 0, the apocentric dis- 
tance is roughly comparable to the virial radius of the host 
system, and that the circular velocity of the host does not 
change significantly with redshift, we estimate that halos 
accreted at z = 10, 5, 1 have been able to complete 83, 32, 
and 5 orbits, respectively, in a ACDM cosmogony. 

Considering that it takes about 10 orbits for a satel- 
lite to lose ~ 95% of its original mass in an orbit with 
fper = 1-5 R s , we conclude that it is unlikely that systems 
accreted after z ~ 2.1 have experienced this much mass 
loss if placed on a comparable orbit. For orbits that ven- 
ture closer to the center than ~ 1.5 R s (<~ 15 kpc for a 
galaxy-sized halo, ~ 0.15 Mpc for a galaxy cluster), halos 
must have been accreted even later in order to retain ;> 5% 
of their initial mass by the present day. 

This calculation emphasizes that substructure is contin- 
ually evolving, and that the surviving population of sub- 
structure halos may be significantly biased towards satel- 
lites with large pericentric distances, which can survive 
for many orbital times, and towards recently accreted ha- 
los, which have yet to be experience significant mass loss. 
A successful model for substructure thus require a good 
understanding of the statistics of accretion events, of the 
orbital parameter distribution of accreting halos, as well 
as a reliable model for their survival time. We explore be- 



low whether simple models of tidal stripping can be used 
to estimate accurately the effects of tides on the evolution 
of the bound mass of substructure halos. 

4. MODELING MASS LOSS 

Over the years, a number of sophisticated theoretical 
prescriptions for tidally-driven mass loss have been devel- 
oped, largely within the context of the evolution of globu- 
lar clusters in external tidal fields (e.g., Gnedin et al. 1999 
and references therein; see also Taylor & Babul 2001). 
However, these techniques (e.g., Fokker-Planck models) 
are computationally expensive to implement and not par- 
ticularly well suited to the analysis of the evolution of sub- 
structure in collisionless systems. Because of this, most 
work on the evolution of substructure has used relatively 
simple prescriptions based on the tidal approximation and 
on the impulse approximation to predict the effects of tides 
on subhalos (Moore et al. 1996, Ghigna et al. 1998, 2000, 
Tormen et al. 1998, K99a). We examine in this section the 
applicability of simple models based on these approxima- 
tions with the aim of deriving a simple prescription that 
traces accurately, and as a function of time, the mass loss 
from subhalos driven by external tides. 

4.1. Tidal Radii 

The simplest approach to mass loss modeling is based 
on the 'tidal approximation,' which assumes that all mass 
beyond a suitably defined tidal radius is lost in a single or- 
bit. The tidal radius, r t , is usually defined as the distance 
from the center of the satellite beyond which the differen- 
tial tidal forces of the host potential exceed the self-gravity 
of the satellite. The precise definition of the tidal radius 
depends on a number of assumptions. For example, if we 
assume that the gravitational potential of both the satel- 
lite and the host are given by point masses, and that the 
satellite is small compared to its distance from the center 
of the host, we find the definition of tidal radius familiar 
from the Roche limit, 



/ m \ 1 / 3 



(4) 



where m is the mass of the satellite, M is the mass of 
the host, and R is their relative distance. A more general 
expression may be obtained by considering the extended 
mass profiles of both systems (see, e.g., Tormen et al. 1998, 
K99a), 
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A further complication is introduced by the centrifugal 
force experienced by the satellite. If this effect is consid- 
ered for a satellite on a circular orbit we find, for point 
masses, the Jacobi limit, 
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According to eqs. 4 and 6 the mean density of the satellite 
within the tidal radius is three (two) times the mean den- 
sity of the host within R in the Jacobi (Roche) limit, and 
less than twice the mean density of the host in the more 
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general case of eq. 5. 5 Although these formulae assume 
circular orbits, they are usually applied to eccentric orbits 
assuming that R corresponds to the pericentric distance, 
where tides are strongest. Unless otherwise specified, we 
will refer to the tidal limit, r t , as the radius defined by 
cq. 5. Wc prefer this definition because it takes into ac- 
count the extended mass profiles of the satellite and host, 
and also because it has been used in previous analyses of 
substructure (Tormen et al. 1998, K99a) and therefore it 
is interesting to test its accuracy. 

Figure 6 shows the tidal radius (left panel), as well as 
the satellite mass within the tidal radius (right panel), as 
a function of orbital radius for a satellite-host system with 
M^ost = 300 m sat . Three different estimates of the tidal 
limit are shown, the smallest being the Jacobi limit, the 
largest being the radius defined by eq. 5. The various es- 
timates of the tidal radius vary from r t ~ 1 r s to 10 r s for 
orbital radii in the range R ~ 1 — 10 R s . Note that in all 
cases the tidal radius becomes comparable to the minimum 
binding radius, rbi nc i, discussed in § 2, only for orbits that 
venture well inside R — R s , so we would expect from the 
simple argument presented in § 2 that none of the orbits 
with Mhost = 300 m sat listed in Table 2 should lead to full 
disruption. 

4.2. Total disruption of NFW halos 

Do substructure halos subject to strong tides fully dis- 
rupt? In order to test the conditions necessary for total 
disruption, we carried out a number of N = 10 5 simula- 
tions where the tidal radius is comparable to or smaller 
than the minimum binding radius rbi n d ~ 0.77 r s . These 
are the circular orbits listed in Table 2 with orbital radii 
T'circ = 1-5 R s , and Mh os t = 300 to 3000 m sat . The evo- 
lution of the bound mass of these satellites is shown in 
Figure 7 by solid curves labelled by the value of the pa- 
rameter r t /r s . The dashed curves correspond to N = 3000 
satellites with tidal radii significantly larger than rbind- 
Figure 7 shows clearly that, for r ( < 2r s ~ 2.6rb; n d the 
satellite halo fully disrupts in just a few orbits. Satellites 
with larger tidal radii, on the other hand, lose mass at an 
ever decreasing rate, and may actually survive indefinitely 
the effects of tides. 

This result appears to contradict the findings of Moore 
et al. (2001) who claim that "the singular cores of sub- 
structure halos always survive complete tidal disruption 
although mass loss is continuous and rapid.' ' This conclu- 
sion was based on simulations of a Hernquist (1993) model 
halo on a circular orbit within a Galactic potential, with 
r t = a, where a is the scale radius of the Hernquist model, 
where the circular velocity peaks. The minimum bind- 
ing radius of a Hernquist model occurs at rbind = 0.4 a. 
Therefore, the most extreme orbit explored by Moore et 
al (2001) has r t w 2.5 rbind, comparable to the r t /r s = 2.1 
run shown in Figure 7. In the latter case, a small core 
of 0.2% of the original mass remains self-bound even after 
10 orbits. Moore et al find that 0.3% of the original mass 
remains bound after about 5 orbits, consistent with their 
slightly smaller tidal radius. As shown in Figure 7, smaller 
tidal radii lead to rapid and full disruption. We conclude 
that total disruption of NFW satellites does in fact occur, 

5 For satellites with orbital radii R < 2.2R S , K99a suggest a smaller estimate of rt than the one given by eq. 5 in order to account for 
tidal stripping due to resonances between the tidal force of the host and the self-gravity of the satellite. Of the orbits we have simulated, this 
correction applies only those with r per or r c ; rc equal to 1.5 R s and is even then only a minor (7%) correction. 



albeit only for tidal radii smaller than ~ 2 rbind : Moore et 
al concluded that halos survive indefinitely only because 
their simulations did not probe tidal radii small enough to 
ensure disruption. 

4.3. Mass Loss and Tidal Approximation 

How well does the tidal approximation predict mass loss 
during a single orbit? This is shown in the left panel of Fig- 
ure 8, where we compare the mass fraction lost after each 
pericentric passage with that beyond r t , computed using 
the structure of the satellite measured at the preceding 
apoccnter. The tidal approximation is in reasonably good 
agreement with the simulations in the case of strong mass 
loss (;> 10%), but significantly underestimates the effect of 
tides in the case of more moderate mass loss. Comparing 
the open and filled symbols in Figure 8 shows that this 
conclusion is insensitive to the number of particles used in 
the simulations. 

The tidal approximation also significantly underesti- 
mates the effects of tides for circular orbits. This is shown 
in the right panel of Figure 6, where the symbols connected 
by dotted lines indicate the bound mass after completing 
1, 5, and 10 circular orbits at various radii. Even the 
smallest estimate of the tidal radius, i.e. the Jacobi limit, 
underestimates the mass lost in the first orbit by ~ 25% for 
satellites at R < 10 R s . After 10 orbits, satellites placed 
at R — 10 R s have lost <~ 65% of their mass, although 
the tidal approximation would predict that no mass loss 
should occur. We conclude that adopting the tidal ap- 
proximation to predict the effects of tides on substructure 
halos may lead to substantial underestimation of the mass 
loss actually incurred by these systems. 

4.4. Impulse Approximation 

An alternative approach to estimating the effects of tidal 
mass loss is to use the impulse approximation to compute 
the changes in velocity expected for each particle during 
a pericentric passage. The procedure assumes that tides 
operate on a timescale short compared with the internal 
crossing time of the satellite, and is therefore best suited 
for highly eccentric orbits, where a satellite spends a short 
fraction of its orbital time near pcriccntcr. 

We follow a similar procedure to that outlined by 
Aguilar & White (1985). Freezing the configuration of 
a satellite at apocenter, we compute for each particle the 
total change in velocity expected after one orbit if it were 
to move on a trajectory similar to that followed by the 
center of mass of the satellite, 



Av, 



-/ ■ 
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ai(t)dt, 



(7) 



where (t) is the acceleration on particle i due to the host 
potential at each point along the orbit of the satellite. Af- 
ter modifying each particle's velocity in this manner, we 
compute the stripped mass by applying the same proce- 
dure described in § 3.2. 

Unlike the tidal approximation, this method incorpo- 
rates both spatial and kinematic information about the 
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halo particles. As seen in the right panel of Figure 8, this 
extra information results in better agreement between the 
mass loss predicted per orbit and the results of our numer- 
ical experiments. The accuracy of the estimate improves 
as the eccentricity of the orbit increases, as expected since 
the assumptions of the impulse approximation are better 
justified in this case. 

4.5. Cumulative Mass Loss 

The success of the impulse approximation and, to a 
lesser extent, of the tidal approximation, in predicting ac- 
curately the mass stripped in a single orbit relies heavily 
on knowledge of the structure of the satellite before peri- 
centric approach. The results shown in Figure 8 use the 
actual structure of the satellite at each apocenter to pre- 
dict the tidal mass loss during the subsequent pericentric 
passage. Without this knowledge, the accuracy of both 
analytic techniques is severely curtailed. 

This can be seen in Figure 9, where we compare the re- 
sults of applying the impulse approximation with two dif- 
ferent assumptions: (i) using dynamical information from 
the actual simulation at each apocenter (filled circles), and 
(ii) repeatedly applying the impulse approximation to the 
original (unevolvcd) equilibrium configuration of the satel- 
lite (step-like solid lines without symbols). The results 
of the numerical simulations are shown as dotted lines. 
Clearly, unless the structural evolution after each episode 
of mass loss is taken into account, the applicability of this 
analytic technique to the study of the mass loss over sev- 
eral orbits is quite limited. The reason for this is that 
the structure of the bound system evolves away gradually 
from the original NFW profile as tides operate; we turn 
our attention to this issue next. 

5. THE INTERNAL STRUCTURE OF STRIPPED NFW 
HALOS 

The structure of the bound remnant of the satellite 
changes steadily as the satellite gradually loses mass 
through tides. This is shown in Figure 10, where we 
plot the density, mass, and circular velocity profiles of the 
bound remnant in simulations with 10 5 particles. These 
profiles are shown at apocenter, when the satellite is clos- 
est to equilibrium, and illustrate a number of interesting 
trends. 

Firstly, as expected, much of the tidally stripped mass 
is lost from the outer regions of the halo, steepening the 
outer density profile. The effects of tidal heating are, how- 
ever, not restricted to the outer regions. As satellites lose 
mass, the density in regions closer to the center also de- 
creases significantly. Note that tides do not apparently 
lead to the formation of a constant-density core near the 
center. These effects tend to reduce the peak circular ve- 
locity, Vr„ ax , of stripped satellites and to shift the location 
of the peak, r max , inwards. The correlation between the 
peak circular velocity of stripped halos and their remain- 
ing bound mass is shown in the top panel of Figure 12. 

1/3 

We find that V m&K scales roughly as m h ' nd as indicated by 
the solid line in that figure. 

The second trend to note is that the structure of a 
stripped halo seems to be dictated by a single parame- 
ter: the total amount of mass lost. This is clear from the 
mass- and velocity-profile panels in Figure 10, which show 
a gradual and smooth progression in the radial distribu- 



tion of mass as a function of remaining self-bound mass. 
Note that this trend is relatively independent of the or- 
bital parameters of the satellite, i.e., the mass profile of a 
satellite that has lost almost 60% of its initial mass over 
five orbits is much the same as that of a satellite that has 
lost the same mass fraction in a single orbit by passing 
much closer to the centre of the host potential. 

It is possible to describe the structure of a stripped halo 
through a simple modification of the NFW profile (eq. 1), 

ft 

P(r) = 1 + ( r /r te )3 PNFW (r) ' (8) 

where ft is a dimensionless measure of the reduction in 
central density, and r te is an 'effective' tidal radius that 
describes the outer cutoff imposed by tides. For f t = 1 
and r te ^> r s , eq. 8 reduces to the original NFW formula. 

Fits to various mass profiles of stripped halos using eq. 8 
are shown in Figure 1 1 and are seen to describe accurately 
the radial mass structure of these systems. The parame- 
ters ft and rte are not independent, and are both deter- 
mined by the mass fraction of the satellite that remains 
bound, m bnd , as shown in Figure 12. The effective tidal 
radius may be computed from the bound mass using the 
polynomial fit, 

logr te = 1.02 + 1.38 logm bnd + 0.37 (logm bnd ) 2 , (9) 

whereas ft follows directly from r te by imposing the con- 
dition that the total mass be mbnd- In practice, ft is well 
approximated by 

log f t = -0.007 + 0.35 log m bnd 

+ 0.39 (logm bnd ) 2 + 0.23 (logm bnd ) 3 , 

as can be seen from the fits shown in Figure 12. 

Finally, the density-profile panel in Figure 10 suggests 
that, after a few orbits, tides tend to impose a well de- 
fined outer cutoff in the mass distribution. This can also 
be seen in Figure 12, where we see that the simulations 
with r pcl = 6.0 R s and 3.0 R s tend asymptotically to a well 
defined minimum r te (and bound mass, see Figure 9). The 
fper = 1-5 R s run, on the other hand, does not participate 
in this trend, as it appears to approach asymptotically 
total disruption. 

5.1. Cumulative Mass Loss Revisited 

We can use eq. 8 to predict fairly accurately the struc- 
ture of the satellite after each episode of mass loss, and 
use this to address analytically the cumulative mass loss 
history of a satellite halo. As discussed in § 4.3, this im- 
plies that we can apply either the tidal or the impulse 
approximation one orbit at a time in order to predict the 
tidal mass loss over several orbital periods. Our prescrip- 
tion for predicting cumulative mass loss using the impulse 
approximation is as follows: 

• apply the impulse approximation to an equilibrium 
realization of an NFW halo and compute the mass 
lost after one orbit using the unbinding procedure 
described in § 3.2, 
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• use the remaining bound mass to compute the mod- 
ified halo profile parameters f t and r te needed to 
calculate the mass profile of the stripped satellite 
(cqs. 8, 9 and 10), 

• construct an isotropic, equilibrium halo model with 
this new mass profile, 

• repeat the whole procedure using the new stripped 
halo model 

In the case of the tidal approximation, the procedure is 
more efficient, since we calculate the tidal radius using the 
analytic NFW and stripped halo mass profiles instead of 
generating N-body realizations of halo models after each 
orbit. 

The success of these methods may be judged from Fig- 
ure 9, where we show the predicted evolution of the bound 
mass for the three orbits shown. The impulse approxima- 
tion (open circles) reproduces the results of the numerical 
simulations reasonably well, especially in the first 4 — 5 
orbits. The predictions deviate from the simulation re- 
sults at later times, primarily due to discrepancies in the 
precise structure of the stripped halos; the change in the 
shape of the profile caused by tidal heating is not perfectly 
captured by the stripped halo profile of eq. 8. A similar 
result is obtained by applying the tidal approximation to 
compute the fractional mass loss during each orbit (open 
squares) . 

In conclusion, combining the evolving structure of 
stripped halos (eq. 8) with simple estimates of the mass 
loss per orbit provides a simple theoretical tool that may 
be used to guide the analysis and interpretation of sub- 
structure studies in cosmological simulations. 

6. SUBSTRUCTURE IN THE MILKY WAY 

The results from the previous sections indicate that sub- 
structure halos are fully disrupted only if they venture 
close to the center of their host systems and, even then, 
disruption takes a number of orbits to complete. It is thus 
not surprising that high-resolution cosmological N-body 
simulations have been able to resolve a wealth of sub- 
structure within virializcd dark matter halos made up of 
the partially disrupted remains of earlier accretion events 
(Moore, Katz & Lake 1996, K99a). 

What are the luminous counterparts of substructure ha- 
los? In a galaxy cluster, it is tempting to identify sub- 
structure halos with dynamical tracers of the individual 
galaxy population, but in the case of a galaxy-sized sys- 
tem such association would imply that several hundred 
satellites should orbit the Milky Way. This is in stark 
contrast with the dozen or so known Milky Way satellites, 
and implies that most subhalos must have failed to form 
a significant number of stars if these models are to match 
observations (Kauffmann, White & Guiderdoni 1993, M99, 
K99b). 

We illustrate this in Figure 13, where we plot, with 
dashed lines, the (peak) circular velocity function of sub- 
structure halos in three high-resolution simulations of 
galaxy-sized CDM halos, as compiled by Font et al. (2001). 
Scaled to the virial velocity of the host halo, V200 = 
K(f20o)j this function is roughly independent of the mass 
of the host and of the value of the cosmological param- 
eters, making it a robust theoretical prediction ideal for 



comparison with observations of Galactic satellites (K99b, 
M99, Font et al. 2001). 

Such comparison requires, however, good estimates for 
quantities that are not directly accessible to observation: 

(i) the virial velocity, V200, of the Milky Way's halo, and 

(ii) the peak circular velocity of halos surrounding the 
Galactic satellites. Conclusions from this exercise are 
thus sensitive to assumptions made to infer such quan- 
tities from observables. For example, one may assume (as 
in M99) that the Milky Way's halo has a virial velocity 
of V200 = 220 km s -1 , and that peak halo circular veloc- 
ities in dwarf irregular (dlrrs) satellites are identical to 
the rotation velocity of their gaseous disks. In addition, 
peak circular velocities are inferred for dwarf spheroidals 
(dSphs) by assuming that stars in these systems are on 
isotropic orbits in isothermal potentials. 

The result of these assumptions is shown in Figure 13 
by the filled circles joined with a solid line. The figure in- 
cludes only satellites within the virial radius of the Milky 
Way, r 2 oo ~ (V2oo/km s _1 )/i _1 kpc <~ 338 kpc. With 
these assumptions, M99 concluded that dSphs populate 
halos with peak circular velocities of order ~ 10 kms -1 
(V max ~ 0.05 V2oo)- According to Figure 13, one expects 
several hundred substructure halos of comparable veloc- 
ity, implying that fewer than one in ten of such halos are 
inhabited by luminous dSphs. 

What makes the very few systems that host luminous 
dSphs distinct from the rest? This question has prompted 
suggestions that dwarf spheroidals formed only in sub- 
halos that collapsed before the universe was fully reion- 
ized, since late-collapsing systems would experience dif- 
ficulty retaining and cooling their baryonic component 
after reionization (Bullock, Kravtsov & Weinberg 2000, 
Somerville 2001, Benson et al. 2001). 

White (2000), on the other hand, has questioned this 
scenario on the grounds that the isothermal assumption 
used by M99 to derive peak circular velocities for dSphs 
is not supported by the results of direct numerical simu- 
lations. As discussed in § 2, numerical simulations indi- 
cate that halo circular velocities decrease systematically 
towards the center, implying that if stars populate the in- 
nermost regions of subhalos, then the isothermal assump- 
tion may substantially underestimate the peak circular ve- 
locities of their surrounding halos. Assuming that dSphs 
are surrounded by NFW halos, White (2000) estimated 
that dSphs may plausibly inhabit potential wells with cir- 
cular velocities up to a factor of 3 times larger than in- 
ferred under the isothermal assumption. A correction of 
this magnitude would reconcile, at the high mass end, the 
Milky Way satellite velocity function with the subhalo ve- 
locity function (shown with dashed lines in Figure 13). 
We emphasize, however, that the magnitude of the correc- 
tion depends sensitively on the inner structure of subhalos, 
which is poorly resolved even in the best cosmological sim- 
ulations available at present. 

Resolving this controversy is important for a couple of 
reasons. Firstly, if dSphs indeed inhabit halos with peak 
circular velocities as high as SOkms" 1 (~ 15% of V200), 
then essentially all such halos would host a dSph and the 
scarcity of dSphs would just reflect the relative scarcity of 
massive substructure halos. This would suggest that reion- 
ization has not played a major role in hindering the for- 
mation of dwarf spheroidal galaxies. The second reason is 
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that if dSphs are indeed surrounded by halos that massive, 
it would be easier to understand how they have managed 
to retain their gaseous components so as to undergo several 
distinct episodes of star formation (Mateo 1998). We use 
below our results for the structure of stripped substructure 
halos in order to place constraints on the peak circular ve- 
locity of dark halos surrounding dwarf spheroidals. 

6.1. The Dark Halos of Dwarf Spheroidals 

If the shape of the dark halo mass profile is specified, it is 
possible to use the velocity dispersion and spatial distribu- 
tion of the stellar component to determine the structural 
parameters of halos surrounding dSphs. Assuming spher- 
ical symmetry and that the mass profile of substructure 
halos can be approximated by eq. 8, we may use Jeans' 
equations to estimate the halo parameters. The velocity 
dispersion, 0% , of stars embedded in a system with circular 
velocity profile V c (r) is given by 

where M*(r) is the stellar mass profile. As is customary, 
we use theoretical King models (King 1966) to describe 
the dSph luminosity profiles; the parameters used for each 
system are listed in Table 3, together with the appropriate 
references. The analysis assumes a stellar mass-to-light ra- 
tio of 2 (solar units) in the F-band, although we note that 
our conclusions are rather insensitive to plausible choices 
of this parameter. 

The dark halo contribution to the circular velocity pro- 
file in cq. 11 is fully specified once three quantities are 
defined: the physical parameters of unstripped NFW ha- 
los, p s and r s (see eq. 1), as well as the total fraction of 
mass lost (if any). The NFW parameters depend on cos- 
mogony, and have been selected assuming a ACDM cos- 
mogony (fl = 0.3, Oa = 0.7, h = 0.65, og = 0.9), follow- 
ing the prescription of Eke, Navarro & Stcinmetz (2001). 

6.1.1. Carina 

The solid curves in the top left panel of Figure 14 il- 
lustrate the circular velocity profile of (stripped) halos 
chosen to match the observed stellar velocity dispersion 
(cr* w 6.8km s _1 ) of the Carina dwarf spheroidal. Different 
curves correspond to various assumptions about the frac- 
tion of mass that has been lost to stripping, ranging from 
95% to a case where no mass loss has occurred (thick solid 
line) . The dashed curves in the same figure show the circu- 
lar velocity profiles of the corresponding halos before strip- 
ping. The NFW concentration parameter (c = r 2 oo/?"s) of 
the halos before stripping varies from c ~ 14 in the case 
of no mass loss to c ~ 12 in the case of 95% mass loss. An 
upward vertical arrow labelled r t i indicates the location of 
the 'tidal' cutoff in the luminous profile, obtained from the 
literature (see Table 3). Downward vertical arrows indi- 
cate the location of the tidal radii, r te , corresponding to 
each stripped halo. 

Figure 14 shows that the stellar velocity dispersion con- 
strains tightly the circular velocity within the 'tidal' cutoff 
of the stellar component, V c (ru), regardless of the degree 
of stripping. Indeed, all circular velocity curves that match 
this constraint (solid lines in Figure 14) agree closely at r t i , 
regardless of the stripped mass fraction, and have roughly 



the same peak circular velocity. This implies that the peak 
velocity of Carina's halo is well defined within this model, 
and is unlikely to be less than 54kms _1 , well in excess of 
its stellar velocity dispersion of 6.8 kms -1 . 

An important consequence of this conclusion is that the 
true tidal radius of Carina's halo should greatly exceed the 
cutoff in the surface brightness profile that is usually in- 
terpreted as a tidal limit. Even assuming that 95% of the 
original halo mass has been stripped the true tidal radius is 
still ~ 20 times larger than derived from King model fits to 
the luminosity profile. This suggests that the stars beyond 
the luminous cutoff detected, for example, by Majewski et 
al. (2000), may not actually represent a true population of 
unbound, extra-tidal stars but rather correspond to a radi- 
ally extended component bound to Carina. Confirmation 
of the extra-tidal nature of such stars through independent 
means, such as detecting tidal tails, or apparent 'rotation' 
in Carina's extended envelope (see Majewski et al. 2000 
for details) would thus provide a strong argument against 
this conclusion. 

We note that Mayer et al (2001) have recently reached a 
similar conclusion, noting that if Carina's luminous cutoff 
is truly tidal in origin this would be very difficult to rec- 
oncile with the massive dark halos expected to surround 
dwarf spheroidals. Mayer et al's conclusion relies on their 
identification of dSphs as tidally-stirred dlrrs; our results, 
on the other hand, are independent of the formation mech- 
anism for dSphs and rely solely on the assumption that 
stars are in dynamical equilibrium within a dominant dark 
matter halo, an assumption well founded on current inter- 
pretation of the observational evidence. Our conclusion is 
also supported by the work of Stoehr et al (2002), who ar- 
gue that the internal structure and kinematics of the Milky 
Way satellites are in excellent agreement with substructure 
found in high-resolution CDM simulations provided that 
the satellites inhabit the most massive substructure halos. 

6.1.2. Draco 

The above conclusion may also apply to other dSphs, 
where apparent radial cutoffs in the outer surface bright- 
ness profiles could actually be of little dynamical signifi- 
cance and may not necessarily correspond to actual tidal 
features. One especially interesting case is Draco. Us- 
ing the core and tidal radii derived by Irwin & Hatzidim- 
itriou (1995, hereafter IH95) from photographic mate- 
rial (r coro = 158 pc, ru — 498 pc), our analysis re- 
quires a very large halo peak velocity, Kn ax ~ 130 km s -1 , 
to accommodate the observed stellar velocity dispersion 
(cr* ~ 9.5 kms -1 , Mateo 1998). This uncomfortably large 
velocity, almost comparable to that typically associated 
with L* galaxies, is difficult to understand given the low 
luminosity of a dwarf system such as Draco. 

One possibility is that gas dynamics and star forma- 
tion have significantly altered the mass profile of the 
halo so that it is no longer well-described by an NFW 
model. Assuming an isothermal halo model with a flat 
circular velocity profile as in M99, for example, yields 
Knax = v^er* = 13.4 km s -1 for Draco. Given the as- 
sumption of NFW halo models, however, the main reason 
for the large peak velocity is the rather small tidal radius 
derived for Draco by IH95; this implies a very large dark 
matter density within r t i to account for the observed ve- 
locity dispersion, which in turn requires a very massive 
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halo. (Within a given radius, more massive NFW halos 
are also denser, as shown by the dashed lines in Figure 14; 
see also NFW and Eke, Navarro & Steinmetz 2001.) 

A similar conclusion was reached by Burkcrt (1997), who 
argued that the mean density within r t i in Draco is much 
higher than the mean Galactic density interior to Draco's 
distance to the center of the Galaxy and that, therefore, 
the cutoff radius derived by IH95 cannot represent a true 
limitation imposed by tides. We emphasize, however, that 
our conclusion is independent of Burkert's analysis, and 
is based solely on our results for the structure of stripped 
NFW halos, as discussed above. 

Our results thus support the revision to Draco's lumi- 
nosity profile derived by Odcnkirchcn ct al. (2001) on the 
basis of data from the Sloan Digital Sky Survey (SDSS). 
These authors fail to find a well defined tidal feature in 
Draco's surface brightness profile and conclude, by fitting 
empirical King profiles (King 1962) that, if one exists, it 
must lie beyond 40.1 arcmin, more than 40% farther from 
the center than IH95's original determination of 28.3 ar- 
cmin (along the major axis). Fitting theoretical King mod- 
els (King 1966) results in an even larger limiting radius, 
fti J> 49.4 arcmin. 

Repeating our analysis with the SDSS theoretical King 
model parameters (r coro = 179 pc, r t i = 1020 pc and as- 
suming the distance to Draco as Odenkirchen et al. (2001) 
(d = 71 ± 7 kpc), we find a peak circular velocity of or- 
der 79 kms -1 for Draco, more in line with what might be 
expected from the very low luminosity of this system (top 
right panel of Figure 14). The concentration parameter of 
the Draco model halos varies from c ~ 14 in the case of no 
mass loss to c ~ 11 in the case of 95% mass loss. Interest- 
ingly, we also find good agreement with the recent work of 
Kleyna et al. (2001), who find a mean mass-to-light ratio 
within ~ 30 arcmin of (330 ± 125)M /L Q . Within the 
same radius, our model gives M/L — 346 ± 128M /L Q . 

6.2. Dependence on Accretion Redshift 

As the discussion above illustrates, the luminosity pro- 
file and velocity dispersion of dSphs constrain tightly the 
circular velocity within the luminous boundary of the sys- 
tem, V c (r t i). For unstripped NFW halos identified at a 
given redshift, this constraint is sufficient to determine the 
virial mass of the halo: our analysis in § 6.1 assumes NFW 
parameters appropriate for halos identified at z = and 
shows that the V c (ru) constraint can only be matched if 
fairly massive halos surround dSphs. Since tides remove 
mass from all radii, the mass of unstripped NFW halos 
that match V c (r t i) provide a lower limit to the mass of 
halos where dSphs formed. 

It is possible, however, that the characteristic densi- 
ties and scale radii of dSph halos are better described by 
the parameters of NFW halos identified at higher redshift. 
This is because substructure halos are expected to evolve 
differently from the relatively isolated systems dealt with 
by the numerical work on which our estimates of p s and 
r s are based (see, e.g., NFW, Eke, Navarro & Steinmetz 
2001). In particular, subhalos effectively stop accreting 
mass after being incorporated into a more massive sys- 
tem, implying that their structural parameters may very 
well reflect the NFW concentrations prevalent at the red- 
shift of accretion rather than at present. 

How does this affect the peak circular velocity derived 



for dSph halos? Dark halos are denser at higher redshift, 
so it is in principle possible for lower mass halos to match 
the V c (r t ) constraint. This is shown for Carina in the bot- 
tom left panel of Figure 14, which is identical to the top 
left panel but adopting NFW parameters appropriate for 
halos identified at z = 5 rather than at z = 0. The con- 
centration parameter is c <~ 3.0 both cases. The effect is 
noticeable, but weak; the minimum peak velocity compat- 
ible with the V c {r t i) constraint for Carina is 53.8 kms^ 1 
for z = 0, and 41.6 kms^ 1 for z — 5. The bottom right 
panel of Figure 14 shows the same effect for Draco. The 
concentration parameter here varies from 2.9 in the case of 
no mass loss to 2.5 in the case of 95% mass loss. The min- 
imum peak velocity is 79 kms -1 for z = 0, and 60 kms -1 
for z = 5. 

Thus, a firm lower limit on the peak velocity of dSph 
halos requires an upper limit on the accretion redshift of 
dSphs into the Milky Way's halo. The accretion redshift 
itself is highly uncertain, but it is possible to derive a rough 
upper limit by considering the present-day Galactocentric 
distance of a dSph. As discussed in § 3.2, the spherical 
infall model suggests that satellites accreted early should 
orbit closer to the Galaxy than late accreting systems. In 
other words, the present-day Galactocentric distance of a 
satellite provides a rough upper limit to the virial radius 
of the Galaxy at the time of accretion. 

Given the one-to-one correspondence between virial ra- 
dius and accretion time implied by the spherical secondary 
infall model (§ 3.2), one can use this assumption to derive 
accretion redshifts for all dSphs. For example, a satellite 
presently at the virial radius would be assigned an ac- 
cretion redshift of z acc — 0; the same procedure assigns 
Zacc ~ 3.3 to Ursa Minor, the closest dSph at 68 kpc, and 
z acc ~ 0.6 to Leo I, the farthest in our sample at 254 kpc. 
These are best regarded as upper limits to the true accre- 
tion redshift: the true apocenter may lie well beyond the 
present-day Galactocentric distance. Likewise, the semi- 
major axis of the satellite's orbit may have been eroded by 
dynamical friction. Both effects would reduce the actual 
accretion redshift relative to our estimate. 

6.3. The Circular Velocity Function of dSph Halos 

We list in Table 3 the upper limits to the accretion red- 
shift (z acc ) derived for each dSphs in our sample using 
the procedure described in the previous subsection. These 
can in turn be used to derive firm lower limits to the 
peak circular velocity of their surrounding halos, listed 
as V max (z acc ) in Table 3. To be conservative, we adopt 
z aC c=10 for all dSphs (well in excess of the upper limits 
shown in Table 3) and derive the circular velocity func- 
tion of Milky Way satellites. This is shown in Figure 13 
as the leftmost boundary of the shaded region; the right- 
most boundary corresponds to assuming z acc = for all 
satellites. 

Figure 13 thus confirms the suggestion of White (2000) 
that in all likelihood the circular velocity of dSph halos 
peaks at values that greatly exceeds the circular speed 
at the luminous boundary V c (r t i), shown as open circles 
connected by a dotted line. This implies that their dark 
halos must extend well beyond r t i- The shaded region is 
intended to illustrate the uncertainty in this function as- 
sociated solely with the accretion redshift dependence of 
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our Vinax estimates (true uncertainties could be substan- 
tially larger). Since these estimates neglect any possible 
stripping, they truly represent lower limits to the actual 
peak circular velocity of dSph halos. 

Even under such conservative assumptions, we can see 
from Figure 13 that, within the uncertainties, there ap- 
pears to be no major discrepancy between the number of 
massive satellites expected in the CDM scenario and the 
known satellite companions of the Milky Way. Only satel- 
lites with peak circular velocities <; 35 km/s (<; 0.16 V200) 
seem to lack readily detectable luminous components. The 
scarcity of satellites around the Milky Way reflects the 
small number of substructure halos with circular veloci- 
ties exceeding ~ 35 km/s, and a simple scenario where 
feedback prevents galaxies from forming in systems below 
a critical circular velocity seems consistent with observa- 
tion. Our results thus indicate that it may be possible to 
explain the abundance of dSphs without invoking a sub- 
stantial role for the reionization of the universe. 

Even under such conservative assumptions, we can see 
from Figure 13 that, within the uncertainties, there ap- 
pears to be no major discrepancy between the number of 
massive satellites expected in the CDM scenario and the 
known satellite companions of the Milky Way. Only satel- 
lites with peak circular velocities <; 35km s _1 (<; 0.16 V200) 
seem to lack readily detectable luminous components. The 
scarcity of satellites around the Milky Way thus reflects the 
small number of substructure halos with circular velocities 
exceeding ~ 35 kms -1 . 

7. SUMMARY 

We have used N-body simulations to investigate the 
structural evolution of substructure halos orbiting within 
a massive host system modeled as a static potential. We 
assume that the equilibrium structure of unstripped halos 
can be approximated by the density profile proposed by 
NFW and that dynamical friction, halo-halo encounters, 
and departures from spherical symmetry play a minor role 
in the evolution. 

Our main conclusions may be summarized as follows. 

1. The structure of N-body realizations of equilibrium 
NFW halo models evolve significantly as a result of 
discreteness effects associated with encounters be- 
tween particles. Such effects lead to the formation of 
an isothermal "core" of size ~ r s on the evaporation 
timescale at the scale radius, i C va P (r B )- At later times, 
collisional effects lead to the progressive "collapse" 
of this core, increasing the mean density within 1 r s 
by 50% after t cvap (r s ). The survival of small- N ha- 
los may be affected by this process in cosmological 
simulations; in particular, the evolution of systems 
with N <^ 200 should be carefully monitored in order 
to ensure that discreteness effects do not unduly bias 
the results. 

2. Satellite halos orbiting in the tidal field of a massive 
host lose mass continuously as a result of tidal strip- 
ping, although the mass loss rate slows down signif- 
icantly as the satellite becomes confined within its 
tidal radius. When the tidal confinement approaches 
the minimum binding radius, r^nd, NFW models are 
fully disrupted. Quantitatively, we find that NFW 



halos on orbits where r t < 2r s <~ 2.6 rbmd at peri- 
center are fully disrupted. 

3. Although tides preferentially strip mass from the 
outer regions, tides also cause the halo to expand 
and the central density to decrease after each peri- 
centric passage. As a result, simple models based on 
the tidal-limit approximation underestimate signifi- 
cantly the total tidal mass loss over several orbits. 

4. The equilibrium structure of stripped NFW halos de- 
pends mainly on the fraction of mass lost, and can 
be expressed in terms of a simple correction to the 
original NFW profile (eqs. 8, 9, 10). Models based 
on the impulse or tidal approximation that take into 
account this evolution in structure are found to ac- 
count reasonably well for the mass lost through sev- 
eral pericentric passages. 

Applying these results to the dwarf spheroidal compan- 
ions of the Milky Way, we conclude that their surrounding 
dark matter halos have circular velocity curves that peak 
well beyond the luminous radius at values significantly 
higher than expected from the stellar line-of-sight veloc- 
ity dispersion. The modeling also suggests that the true 
tidal radius of dSphs may lie well beyond the radial cut- 
off observed in surface brightness profiles, suggesting that 
these arc not really tidal in origin but rather are features 
in the luminous profile of little dynamical relevance. Our 
results for Draco, in particular, strongly suggest that the 
tidal radius should be much larger than the 28.3 arcmin de- 
rived by Irwin & Hatzidimitriou (1995), in agreement with 
the recent Sloan survey data reported by Odenkirchen et 
al. (2001). 

Our results imply a close correspondence between the 
most massive substructure halos and the known satellites 
of the Milky Way and suggests that only substructure ha- 
los with peak circular velocities below 35 km s _1 lack 
readily detectable luminous counterparts. The scarcity 
of satellites around the Milky Way just reflects the small 
number of substructure halos with circular velocities ex- 
ceeding ~ 35 kins" 1 . This implies that it may be possi- 
ble to explain the abundance of dSphs without invoking 
a substantial role for the reionization of the universe in 
hindering the formation of dwarf systems. 

Our conclusions also highlight the fact that the proper- 
ties of galaxies inhabiting low-mass cold dark matter halos 
must vary substantially from halo to halo to be consis- 
tent with observations; vastly different galaxies such as 
the SMC and Draco are predicted to live in halos of not 
very dissimilar circular velocities. Explaining the diversity 
of galaxy properties in low-mass halos remains a challenge 
to hierarchical models of galaxy formation that is unlikely 
to be solved until the effects of mass, collapse epoch, effi- 
ciency of feedback, and environment are properly under- 
stood. 
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was modified to generate our N-body halo models. We 
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IBM/SP3 supercomputer at the University of Victoria. 
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Table 1 

Characteristic timescales for an isolated N = 10 6 NFW Halo 



r 


JVb(r) 


^o/^cross 


tcoll /^cross 


tevap ft cross 


0.025 


100 


120 


0.17 


23 


0.05 


500 


190 


0.80 


110 


0.1 


2500 


450 


4.0 


550 


0.2 


10000 


1000 


20.0 


2640 



Note. — Collisional relaxation (t rc iax) and evaporation (i e vap) timescales in the inner regions of an N — 10 6 NFW halo 
model. The timescale for evolution of the inner mass structure observed in the simulations, to, is in reasonable agreement 
with the evaporation timescale, t evap w 136 t ro ia X . 



Table 2 

Parameters of N-body Simulations 



Eccentric Orbits 



N 


M host /m sat 


r ap /R s 




^orb/^cross 


^circ 


n/r s 


3000, 10 5 


300 


10.0 


1.5 


4.8 


0.50 


3.0 


3000, 10 5 


300 


10.0 


3.0 


5.4 


0.75 


4.9 


3000, 10 5 


300 


10.0 


6.0 


6.9 


0.95 


8.8 








Circular Orbits 








N 


M host /m sat 




f circ / ' Rs 


^-orb/^cross 




n/r s 


10 5 


3000 




1.5 


0.19 




0.87 


10 5 


1500 




1.5 


0.65 




1.3 


10 5 


1000 




1.5 


0.79 




1.6 


10 5 


600 




1.5 


1.0 




2.1 


3000, 10 5 


300 




1.5 


1.4 




3.0 


3000 


300 




3.0 


2.9 




4.9 


3000 


300 




6.0 


6.2 




8.8 


3000 


300 




7.5 


8.1 




10.6 


3000 


300 




10.0 


11.5 




13.2 


3000 


300 




30.0 


46.3 




32.2 



Note. — e c i rc = J/J c i IC (E) is the orbital circularity, defined as the ratio between the orbital angular momentum and 
the angular momentum of a circular orbit with the same energy; t OI b / Across is the orbital period in units of the crossing 
time of the satellite halo particles; and r t is the estimate of the tidal radius given by eq. 5 computed at pericenter. 
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Table 3 

Structural and kinematic properties of Satellites of the Milky Way 



Galaxy 


Ltot 




ru 




D 


%acc 






(1O 6 L ) 


(pc) 


(PC) 


kms -1 


(kpc) 




kms -1 


Carina 


0.43 


177 ±28 


581 ± 86 


6.8 ± 1.6 


103 ±5 


2.3 


53.1 


Draco 


0.26 


158 ± 14 


498 ± 47 


9.5 ± 1.6 


82 ±5 


2.8 


73.1 


Fornax 


15.5 


400 ± 43 


2078 ± 177 


10.5 ± 1.5 


140 ±8 


1.6 


33.0 


Leo I 


4.79 


169 ± 19 


645 ± 87 


8.8 ±0.9 


254 ± 30 


0.6 


59.9 


Leo II 


0.58 


162 ±35 


487 ± 60 


6.7 ± 1.1 


208 ± 12 


0.9 


65.3 


Sculptor 


2.15 


101 ±28 


1329 ± 107 


6.6 ±0.7 


79 ±4 


2.9 


25.4 


Sextans 


0.50 


322 ± 42 


3102 ± 1028 


6.6 ±0.7 


89 ±4 


2.6 


17.0 


Ursa Minor 


0.29 


196 ±24 


628 ± 74 


9.3 ± 1.8 


68 ±3 


3.3 


93.2 



Note. — Column (1) lists the name of the satellite galaxy. Columns (2), (5) and (6) list the U-band luminosity, stellar 
central velocity dispersion and galactocentric distance taken from data presented by Mateo (1998). Columns (3) and (4) 
list the King model core radius and tidal radius derived by Irwin & Hatzidimitriou (1995). Columns (7) and (8) list the 
estimated accretion redshift (based on the galactocentric distance) and peak circular velocity of the dark halo surrounding 
each dwarf spheroidal at that redshift. Note that Odenkirchen et al. (2001) derive much larger values of the King model 
core and tidal radii for Draco, r cor c = 179 pc, r t i = 1020 pc. 
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Fig. 1. — Kinetic (K), potential (W) and total (E = K + W) energy profiles for NFW and isothermal distributions truncated at radius r, 
in arbitrary units. Note that the NFW halo (solid lines) has positive total energy if truncated at radii smaller than r bind ~ 0.77 r s , whereas 
the isothermal halo (dashed lines) remains bound at all radii. 
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Fig. 2. — From top to bottom, the mass fraction enclosed within 8.0, 4.0, 2.0, 1.0, 0.5 and 0.25 r s for N = 3000 (solid) and N = 10 5 (dotted) 
isolated NFW halo models as a function of time. 
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Fig. 3.— Top panel: From top to bottom, the mass within 8.0, 4.0, 2.0, 1.0, 0.5 and 0.25 r s for N = 10 6 (solid blue) and N = 3000 (dotted 
red) NFW halo models as a function of time. Solid black lines show the decline in mass within r = 0.5 and 0.25 r s expected from the local 
evaporation timescale of each halo given by eq. 3. Bottom panel: From top to bottom, time evolution of the mass within 0.2,0.1,0.05 and 
0.025 r a for the N = 10 6 halo, normalized to the initial mass within each radius No(r). Solid black lines show fits with exponential decay 
times equal to the evaporation timescale appropriate for each radius. 
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Fig. 4. — Snapshots of satellite halo particles projected onto the orbital plane of the satellite for the (jV, r ap , r per ) = (10 5 , 10.0 R s , 1.5 -R s )- 
The inner circle indicates the pericentric radius and the large circle denotes the apocentric radius of the satellite. Tidally stripped material 
forms elongated tidal tails and spherical shells which extend for many scale radii beyond the remaining bound mass of the satellite halo. 
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Fig. 5. — Bound mass of a satellite halo as a function of time for three different orbits with apocentric and pericentric radii (r ap and r peT , 
respectively) as shown. Results from low and high resolution simulations are shown together for comparison, with number of particles TV and 
softening length e as indicated. Note the good agreement between the two curves despite a difference of a factor of 30 in mass resolution. 
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Fig. 6. — Left: Tidal radius as a function of orbital radius R in units of the satellite scale radius r a for Mi los t = 300m sa t. Dashed 
(dot-dashed) curve shows the tidal radius corresponding to the Roche (Jacobi) limit definition, Psat/Phost = 2(3), where Psat/Phost ' s the 
mean overdensity of the satellite within the tidal radius relative to that of the host within the orbital radius. Solid curve shows tidal radius 
given by equation (5). Right: Mass within tidal radius for the same three definitions of the tidal radius. Symbols indicate the remaining 
bound mass of 3000-particle satellites after completing 1, 5, and 10 circular orbits at orbital radii R = 1.5, 3.0, 6.0, 7.5, 10.0, and 30.0 R a . Note 
that even the Jacobi limit underestimates mass loss after one circular orbit by ~ 25%. 
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Fig. 7. — Bound mass of satellite halos on circular orbits as a function of the number of orbits completed. Each curve is labelled with the 
satellite's tidal radius in units of the scale radius, r s (see Table 2). Solid curves indicate simulations performed with N = 10 , dashed curves 
correspond to satellites with longer orbital periods simulated with 3000 particles. Satellites with tidal radii < 2.6 r^ind ~ 2r s arc totally 
disrupted in a few orbits, whereas satellites with larger tidal radii may survive self-bound for more than 10 orbits and perhaps indefinitely. 
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Fig. 8. — Left: Mass fraction lost between consecutive apocentric passages compared to predictions made using the tidal approximation. 
Results for N = 3000 and N = 10 5 and various orbits are shown, as labelled. Points above the solid line represent mass loss underestimation by 
the tidal approximation; points below the line correspond to overestimates. Right: Same as left panel for impulse approximation predictions. 
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Fig. 9. — Impulse approximation mass loss predictions compared with results for the N = 10 5 simulations (dotted lines) shown in Figure 5. 
Step-like solid curves show predictions made by repeated application of the impulse approximation. Solid circles show independent predictions 
for each orbit based on the structure of the halo at each preceding apocenter, as in Figure 8. Open symbols show predictions of tidal (squares) 
and impulse (circles) approximations using the modified profile given by eq. 8 to recalculate the structure of the stripped halo after each orbit 
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Fig. 10. — Density, mass, and circular velocity profiles of the bound mass of N = 10 5 satellite halos. Profiles are shown at apocenter for 
three different orbits with r ap = 10.0 R B and r per as shown. Density is given in units of m sa t/r|, mass is in units of m sa t, and velocity is in 
units of yGmsat/fj- Thick solid curves correspond to the initial (unstripped) NFW profiles. 
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Fig. 11. — Modified NFW profile fits to selected halo mass profiles from the middle panel of Figure 10. The dashed line shows the original, 
unstrippcd NFW mass profile. The parameters rt e and ft control the truncation radius and the reduction in central density of the stripped 
profile, respectively. 
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Fig. 12. — The peak circular velocity of stripped halos, Vmax (top panel) in units of the initial peak circular velocity, as well as the 
parameters ft (middle panel) and rte (bottom panel), shown as a function of the remaining bound mass fraction m bnd . Solid curves show fits 

1/3 

to the data given by V max oc m, and eqs. 9 and 10 for top, middle and bottom panels, respectively. 
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Fig. 13. — Cumulative (peak) circular velocity function of substructure halos. Velocities are normalized to the circular velocity of the 
parent halo measured at the virial radius. The short-dashed curves correspond to two galaxy-sized halos formed in an SCDM universe (M99). 
Long-dashed curves correspond to a halo formed in the ACDM cosmogony (Font et al. 2001). The dot-dashed curve correspond to results 
reported by K99b. Curves with symbols correspond to the satellites observed around the Milky Way, using different assumptions to compute 
the peak circular velocity of their surrounding halos. Filled circles: Milky Way satellites assuming that stars in dwarf spheroidals are on 
isotropic orbits in isothermal potentials (as in M99) and that V200 = 220km s for the Galaxy. Open circles: circular velocity at the luminous 
cutoff (r t i) assuming that dwarf spheroidals are surrounded by dark halos with NFW profiles. Shaded region is bounded by peak circular 
velocity function of unstripped NFW halos assuming two extreme accretion redshifts, 2 a cc = 10 and 0. See text for details. 
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Fig. 14. — Circular velocity profiles of unstripped and stripped halos chosen to match the surface brightness profile and the stellar velocity 
dispersion of Carina (left panels) and Draco (right panels). Dashed curves show the circular velocity profiles of the corresponding halos before 
stripping. Upward vertical arrows indicate the luminous cutoff radius r ti . The luminous cutoff radius shown for Draco is that of Odenkirchen 
et al. (2001), who find r t ; = 1020 pc. Downward vertical arrows indicate the tidal radius rte corresponding to each of the stripped halos. Top 
and bottom panels show halo profiles adopting Eke, Navarro & Steinmetz (2001) NFW parameters corresponding to halos identified at z = 
and 5, respectively. 



